Source code for hysop.operator.directional.stretching_dir

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


"""
@file stretching_dir.py
Directional stretching frontend (operator generator).
"""
import sympy as sm

from hysop.constants import DirectionLabels, Implementation, StretchingFormulation
from hysop.tools.htypes import check_instance, to_tuple, to_list, first_not_None
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
from hysop.parameters.scalar_parameter import ScalarParameter
from hysop.core.graph.graph import generated
from hysop.fields.continuous_field import Field
from hysop.topology.topology import Topology
from hysop.topology.cartesian_descriptor import CartesianTopologyDescriptors
from hysop.operator.directional.symbolic_dir import DirectionalSymbolic
from hysop.symbolic.relational import Assignment
from hysop.operator.directional.directional import DirectionalOperatorFrontend


[docs] class DirectionalStretching(DirectionalSymbolic): """ Directional stretching using the symbolic code generation framework. """ @debug def __new__( cls, formulation, velocity, vorticity, variables, dt, C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds, ): base_kwds = first_not_None(base_kwds, {}) return super().__new__( cls, name=name, variables=variables, dt=dt, candidate_input_tensors=None, candidate_output_tensors=None, base_kwds=base_kwds, implementation=implementation, exprs=None, **kwds, ) @debug def __init__( self, formulation, velocity, vorticity, variables, dt, C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds, ): """ Initialize directional stretching operator frontend. Vortex stretching is the lengthening of vortices in three-dimensional fluid flow, associated with a corresponding increase of the component of vorticity in the stretching direction due to the conservation of angular momentum. Solves dW/dt = C * f(U,W) where U is the velocity, W the vorticity and f is one of the following formulation: CONSERVATIVE FORMULATION: f(U,W) = div( U : W ) MIXED GRADIENT FORMULATION: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W where : represents the outer product U:W = (Ui*Wj)ij. . represents the matrix-vector product. ^T is the transposition operator U and W are always three dimensional fields. C is scalar or 3d vector-like of symbolic coefficients. A is a scalar symbolic coefficient. f(U,W) is always directionally splittable. Parameters ---------- formulation: hysop.constants.StretchingFormulation The formulation of this stretching operator: CONSERVATIVE: f(U,W) = div( U : W ) GRAD_UW: f(U,W) = grad(U) . W GRAD_UW_T: f(U,W) = grad(U)^T . W MIXED_GRAD_UW: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W where A is a user given scalar or 3d vector like that contains scalar (floating point) coefficient(s) or parameter(s). The most common mixed gradient formulation is A=0.5. GRAD_UW, GRAD_UW_T and MIXED_GRAD_UW formulations are only valid for incompressible flows, see Notes. velocity: Field Velocity U as read-only input three-dimensional continuous field. vorticity: Field Vorticity W as read-write three-dimensional continuous field. variables: dict Dictionary of fields as keys and topology descriptors as values. dt: ScalarParameter Timestep parameter that will be used for time integration. C: scalar or vector like of 3 symbolic coefficients, optional The stretching leading coefficient C can be scalar or vector like. Contained values should be numerical coefficients, parameters or generic symbolic expressions such that C*f(U,W) is directionally splittable. Here * is the classical elementwise multiplication. Default value is 1. A: scalar symbolic coefficient, optional Should only be given for MIXED_GRAD_UW formulations. ValueError will be raised on other formulations. The linear combination coefficients A is a scalar. Contained value should be a numerical coefficient, a parameter (or a generic symbolic expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W) is directionally splittable. Here * is the classical elementwise multiplication. Default value is 0.5. name: str, optional, defaults to 'stretching'. Name of this stretching operator. implementation: Implementation, optional, defaults to None Target implementation, should be contained in available_implementations(). If None, implementation will be set to default_implementation(). base_kwds: dict, optional, defaults to None Base class keywords arguments. If None, an empty dict will be passed. kwds: Keywords arguments that will be passed towards implementation operator __init__. Notes ----- The MIXED GRADIENT formulation is only valid for incompressible flows (ie. only if div(U)=0). It is obtained by expanding the conservative formulation and simplifying with div(U)=0. The linear combination of gradients comes from this first equation and the fact that (grad(U)-grad(U)^T).W = 0. """ check_instance(formulation, StretchingFormulation) check_instance(velocity, Field) check_instance(vorticity, Field) check_instance(dt, ScalarParameter) check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors) check_instance(base_kwds, dict, keys=str, allow_none=True) check_instance(name, str, allow_none=True) if (velocity.dim != 3) or (velocity.nb_components != 3): msg = "Given velocity is not a 3d vector field." raise ValueError(msg) if (vorticity.dim != 3) or (vorticity.nb_components != 3): msg = "Given vorticity is not a 3d vector field." raise ValueError(msg) name = first_not_None(name, "stretching") base_kwds = first_not_None(base_kwds, {}) C = first_not_None(C, sm.Integer(1)) if isinstance(C, npw.ndarray): assert C.ndim in (0, 1), C.ndim assert C.size in (1, 3), C.size if formulation is StretchingFormulation.MIXED_GRAD_UW: A = first_not_None(A, sm.Rational(1, 2)) if isinstance(A, npw.ndarray): assert A.ndim == 0, A.ndim assert A.size == 1, A.size elif A is not None: msg = f"Formulation is {formulation} but A was specified." raise ValueError(msg) elif formulation is StretchingFormulation.GRAD_UW: A = sm.Integer(1) formulation = StretchingFormulation.MIXED_GRAD_UW elif formulation is StretchingFormulation.GRAD_UW_T: A = sm.Integer(0) formulation = StretchingFormulation.MIXED_GRAD_UW exprs = self._gen_expressions(formulation, velocity, vorticity, C, A) super().__init__( name=name, variables=variables, dt=dt, candidate_input_tensors=( velocity, vorticity, ), candidate_output_tensors=(vorticity,), base_kwds=base_kwds, implementation=implementation, exprs=exprs, **kwds, ) from hysop.operator.base.custom_symbolic_operator import ( CustomSymbolicOperatorBase, ) if not issubclass(self._operator, CustomSymbolicOperatorBase): msg = "Class {} does not inherit from the directional operator interface " msg += "({})." msg = msg.format(self._operator, CustomSymbolicOperatorBase) raise TypeError(msg) def _gen_expressions(self, formulation, velocity, vorticity, C, A): from hysop.symbolic.base import TensorBase from hysop.symbolic.field import div, grad frame = velocity.domain.frame U = velocity.s() W = vorticity.s() lhs = W.diff(frame.time) if formulation is StretchingFormulation.CONSERVATIVE: assert A is None, A UW = npw.outer(U, W).view(TensorBase).freeze() # Freezing is important because sympy will split the derivatives, # ie. makes the stretching term not conservative anymore # (and time integration *WILL* fail) # See Taylor-Green Re=1600 without calling freeze(). rhs = div(UW, frame) elif formulation is StretchingFormulation.MIXED_GRAD_UW: assert A is not None, A rhs = (A * grad(U, frame) + (sm.Integer(1) - A) * grad(U, frame).T).dot(W) else: msg = f"Unknown stretching formulation {formulation}." raise ValueError(msg) rhs = C * rhs exprs = Assignment.assign(lhs, rhs) return exprs
[docs] class StaticDirectionalStretching(DirectionalOperatorFrontend): """ Directional stretching using the symbolic code generation framework. """
[docs] @classmethod def implementations(cls): from hysop.backend.host.python.operator.directional.stretching_dir import ( PythonDirectionalStretching, ) __implementations = {Implementation.PYTHON: PythonDirectionalStretching} return __implementations
[docs] @classmethod def default_implementation(cls): return Implementation.PYTHON
@debug def __new__( cls, formulation, velocity, vorticity, variables, dt, C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds, ): base_kwds = first_not_None(base_kwds, {}) return super().__new__( cls, name=name, variables=variables, dt=dt, formulation=formulation, velocity=velocity, vorticity=vorticity, C=C, A=A, candidate_input_tensors=None, candidate_output_tensors=None, base_kwds=base_kwds, implementation=implementation, **kwds, ) @debug def __init__( self, formulation, velocity, vorticity, variables, dt, C=None, A=None, name=None, implementation=None, base_kwds=None, **kwds, ): """ Initialize directional stretching operator frontend. Vortex stretching is the lengthening of vortices in three-dimensional fluid flow, associated with a corresponding increase of the component of vorticity in the stretching direction due to the conservation of angular momentum. Solves dW/dt = C * f(U,W) where U is the velocity, W the vorticity and f is one of the following formulation: CONSERVATIVE FORMULATION: f(U,W) = div( U : W ) MIXED GRADIENT FORMULATION: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W where : represents the outer product U:W = (Ui*Wj)ij. . represents the matrix-vector product. ^T is the transposition operator U and W are always three dimensional fields. C is scalar or 3d vector-like of symbolic coefficients. A is a scalar symbolic coefficient. f(U,W) is always directionally splittable. Parameters ---------- formulation: hysop.constants.StretchingFormulation The formulation of this stretching operator: CONSERVATIVE: f(U,W) = div( U : W ) GRAD_UW: f(U,W) = grad(U) . W GRAD_UW_T: f(U,W) = grad(U)^T . W MIXED_GRAD_UW: f(U,W) = (A*grad(U) + (1-A)*grad(U)^T) . W where A is a user given scalar or 3d vector like that contains scalar (floating point) coefficient(s) or parameter(s). The most common mixed gradient formulation is A=0.5. GRAD_UW, GRAD_UW_T and MIXED_GRAD_UW formulations are only valid for incompressible flows, see Notes. velocity: Field Velocity U as read-only input three-dimensional continuous field. vorticity: Field Vorticity W as read-write three-dimensional continuous field. variables: dict Dictionary of fields as keys and topology descriptors as values. dt: ScalarParameter Timestep parameter that will be used for time integration. C: scalar or vector like of 3 symbolic coefficients, optional The stretching leading coefficient C can be scalar or vector like. Contained values should be numerical coefficients, parameters or generic symbolic expressions such that C*f(U,W) is directionally splittable. Here * is the classical elementwise multiplication. Default value is 1. A: scalar symbolic coefficient, optional Should only be given for MIXED_GRAD_UW formulations. ValueError will be raised on other formulations. The linear combination coefficients A is a scalar. Contained value should be a numerical coefficient, a parameter (or a generic symbolic expression) such that C*(A*grad(U).W + (1-A)*grad^T(U).W) is directionally splittable. Here * is the classical elementwise multiplication. Default value is 0.5. name: str, optional, defaults to 'stretching'. Name of this stretching operator. implementation: Implementation, optional, defaults to None Target implementation, should be contained in available_implementations(). If None, implementation will be set to default_implementation(). base_kwds: dict, optional, defaults to None Base class keywords arguments. If None, an empty dict will be passed. kwds: Keywords arguments that will be passed towards implementation operator __init__. Notes ----- The MIXED GRADIENT formulation is only valid for incompressible flows (ie. only if div(U)=0). It is obtained by expanding the conservative formulation and simplifying with div(U)=0. The linear combination of gradients comes from this first equation and the fact that (grad(U)-grad(U)^T).W = 0. """ check_instance(formulation, StretchingFormulation) check_instance(velocity, Field) check_instance(vorticity, Field) check_instance(dt, ScalarParameter) check_instance(variables, dict, keys=Field, values=CartesianTopologyDescriptors) check_instance(base_kwds, dict, keys=str, allow_none=True) check_instance(name, str, allow_none=True) if (velocity.dim != 3) or (velocity.nb_components != 3): msg = "Given velocity is not a 3d vector field." raise ValueError(msg) if (vorticity.dim != 3) or (vorticity.nb_components != 3): msg = "Given vorticity is not a 3d vector field." raise ValueError(msg) name = first_not_None(name, "stretching") base_kwds = first_not_None(base_kwds, {}) C = first_not_None(C, sm.Integer(1)) if isinstance(C, npw.ndarray): assert C.ndim in (0, 1), C.ndim assert C.size in (1, 3), C.size if formulation is StretchingFormulation.MIXED_GRAD_UW: A = first_not_None(A, sm.Rational(1, 2)) if isinstance(A, npw.ndarray): assert A.ndim == 0, A.ndim assert A.size == 1, A.size elif A is not None: msg = f"Formulation is {formulation} but A was specified." raise ValueError(msg) elif formulation is StretchingFormulation.GRAD_UW: A = sm.Integer(1) formulation = StretchingFormulation.MIXED_GRAD_UW elif formulation is StretchingFormulation.GRAD_UW_T: A = sm.Integer(0) formulation = StretchingFormulation.MIXED_GRAD_UW super().__init__( name=name, variables=variables, dt=dt, formulation=formulation, velocity=velocity, vorticity=vorticity, C=C, A=A, candidate_input_tensors=( velocity, vorticity, ), candidate_output_tensors=(vorticity,), base_kwds=base_kwds, implementation=implementation, **kwds, )